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ABSTRACT 


A spectral method algorithm is developed for the numerical solution of the full six-dimensional 
Vlasov-Maxwell system of equations. Here, the focus is on the electron distribution function, 
with positive ions providing a constant background. The algorithm consists of a Jacobi 
polynomial-spherical harmonic formulation in velocity space and a trigonometric formulation in 
position space. A transform procedure is used to evaluate nonlinear terms. The algorithm is 
suitable for performing moderate resolution simulations on currently available supercomputers 
for both scientific and engineering applications. 



1. The Relativistic Vlasov-Maxwell Equations 


The dynamics of a high-energy, collisionless plasma are described by the relativistic Vlasov- 
Maxwell equations [1]. These nonlinear equations include special relativistic effects [2] and 
couple the equations of the electromagnetic field (Maxwell’s equations) with the evolution 
equations for single-particle distribution functions (Vlasov equations). In the simplest case, a 
plasma has two species, protons and electrons. Protons have charge e>0, mass m p , and 
distribution function/,,, while electrons have charge -e<0, mass: m e , and distribution function/,,. 
Although non-relativistic Vlasov-Poisson and Vlasov-Maxwell systems have received much 
attention in the distant [3-15] and more recent [16-23] past, relativistic systems appear to be 
somewhat less explored, although there have been linear treatments [24-27], Here, an algorithm 
for a spectral method numerical solution of a fully relativistic, nonlinear Vlasov-Maxwell system 
is developed. The goal is to set up the basic framework necessary for numerical simulation in 
the full six-dimensional case. The purpose of this work is to provide a means to enable the 
stuf ly — and ultimately the prediction — of high-energy charged-particle distributions in space, 
both for intrinsic scientific interest and for optimizing human and robotic space exploration. 


The Vlasov equations for the distribution functions/, (x,p,t) and/.(x,p,f) are 
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To determine the self-consistent fields E and B, the Maxwell equations are needed. 
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The sources present in ( 1 .4) and ( 1 .6) are the charge density p and the current density j: 

(1.7) p = p /} + P e =ei{f p -f e )dp 
(1-8) j = j /; +j, =e j(f p -f e )\dp . 

Additionally, E and B can have external components, as well as self-consistent ones, which can 
serve as external drivers of the coupled Vlasov-Maxwell system. 

The solution of this set of integro-differential equations presents a great challenge, because of 
both the nonlinearity of the couplings and the six-dimensional nature of the solution space of the 
distribution function. Except in special cases, nonlinearity requires the use of computer 
simulation, while the presence of six dimensions has, in the past, pushed such simulations 
beyond the capabilities of then-available computer systems. However, we have now begun to 
move into an era of fast computers with large core memories, and these machines are providing 
the resources necessary to perform simulations of six-dimensional continua with a moderate 
amount of resolution. It is in this context that the following algorithm is presented and we hope 
that this will help further the development of computer simulations of the Vlasov-Maxwell 
system and will lead to a greater understanding of the evolution and distribution of high-energy 
charged particles in space (and other) plasmas. 


2. Nondimensional Form of the Equations 

First, let us nondimensionalize the equations of the previous section. To do this, we denote a 
characteristic length by L , and define 


( 2 . 1 ) 




V' = LV. 


In the previous section. Maxwell’s equations were written in Gaussian form, in which the electric 
field E and magnetic induction B have the same units. To nondimensionalize Maxwell's 
equations, we choose B 0 as a characteristic electromagnetic field strength and n a as a 
characteristic number density, so that the fields and sources are transformed into 
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Here, the dimensional and dimensionless quantities are functions of (x,p,f) and (x ,p J ) , 
respectively. Also, using ( 1 .4) and ( 1 .6), a natural choice for B„ is B 0 = eLn 0 , so we will adopt 
this definition. The nondimensional Maxwell’s equations are 


(2.3) V'-B' = 0 

(2.4) V , -E' = 4tc(p + +P_) 

d B T-,, w 

(2.5) — = -V xE 
dt 

(2.6) ^ = V , xB'-4rc(j + +j_). 


The distribution functions, in turn, take the form: 

(2.7) / + (x , ,p , ,/ / ) = — — - /p( x ^PT), /_ (x', p', t') = ~~~~~ /e (x- PT) • 

n n "o 


The nondimensional Vlasov equation is 


(2.8) ^k + v'.|A ± p ± (r + v'xB').|^ = 0. 
dr ox op 


The constants P+ can be given in terms of L, n 0 , and the classical radius of the electron />: 

_ m e _ e 

(2.9) P_=/i 0 ;>L% p + = P-, r e = j. 

m p m e c 

Unless a different characteristic length L is defined for protons and electrons, the relations in 

(2.9) indicate that the dynamic coupling of the electromagnetic field to the proton distribution 
function /+ is less than that of the field to the electron distribution function/, by a factor of m p /m e 
= 1830. Assuming that there is only one overall characteristic length L, then the dynamic 
evolution off. can be thought of as occurring on a static proton background, whose sole purpose 
is to provide overall charge neutrality, at least for times which are short compared to those 
required for an appreciable evolution of/ + . Here, it is the electron distribution function evolution 
that will be of primary concern and, to this end, we will simplify notation by choosing p. = 1 and 
redefining f. = f. Furthermore, we will henceforth drop all primes and accept that all variables 


3 


occurring in all equations are dimensionless. The nondimensional equations we will be 
concerned with are the following: 


(2.10) — + v • -t- - (E + v x B ) • = 0 

dt dx dp 

(2.11) V • B = 0 

(2.12) V-E = 47tp 

r) R 

(2.13) — = -VxE 

dt 

ri F 

(2.14) — = Vx B -4jt j 
dt 

(2.15) p = \-jfdp 

(2.16) j = — J / \ dp . 


3. Formulation in Velocity Space 

The distribution function depends on p rather than v because the six-dimensional phase space 
volume element dxdp = dxdydzdp x dp^dp : and the distribution function /(x,p) are invariant [2] 
under Lorentz transformations, while dxdv and/(x,v) are not. However, limits on the velocity 
components are - 1 < v k < 1 . while the limits on momentum are - °° < p k < °° . For 
computational purposes, we will work in velocity space, since the associated finite domain is 
more commensurate with the finite numerical structure of a digital computer. 

The dimensionless relation between momentum and velocity is 


To transform from a momentum space formulation to one in velocity space, we need to calculate 
the Jacobian of the transformation, and to apply the chain rule of differential calculus to the 
equations of the previous section. These activities require the following partial derivatives, 
which can be derived from (3.1): 
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Here, d is the dimension of the velocity space, and repeated indices, when used, will denote an 
implicit summation from 1 to d : 

T d 0 

(3.4) v" = VjVj = X v f • 

i=l 


The Jacobian for transforming the integrals in (2. 15) and (2.16) is 


(3.5) 


det 


d Pi 
dv k 



Let us also transform the Vlasov equation (2. 1 0) by using (3.2), along with 


(3.6) g(*,v) =f(x,p). 

Then, equation (2.10) becomes 


(3.7) ^ + v • - ( E + u X 6) ( - {b ik - VjY k )(\~v 2 Y 2 JTT = 0 - 

dt dx 


(Again, repeated indices denote summation.) 


Since no massive particles have a speed of c or greater, we must have g — 0 for v >1. The 
distribution function g can be written as 

- v 2 )" G(x, v), v 2 <l, n> 0 

0, r 2 > 1 

We will assume that G can be represented by an analytic function in the v k over the domain 
i’ 2 < 1 . Also, using n > 1 in (3.8) will ensure that not only the distribution function, but also its 

first derivative, will be zero at v 2 = 1 . 


(3.8) g(x,v) = 


5 



Using (3.5) and (3.8), the density integral (2.15) becomes 

(3.9) p = l- °jfdp = \- J(l-v 2 )"” T_1 G//(l-v 2 )r/v. 

— OO —1 

Here, H{z) is the step function, which is equal to unity if r > 0. and equal to zero otherwise. The 
current integral (2.16) similarly becomes 

(3.10) j = - J/vrfp = - jv(l-v 2 )"' T_1 G//(l-v 2 )r/v 

— OO — 1 

Using (3.9), equation (3.8) becomes 

(3.11) ^ + v~-(E + uxB),(l-r 2 f 
dt dx 

The choice of n in (3.9), (3. 10) and (3. 1 1) is determined so as to simplify the algorithm. 


(S/A- - v i v k ) 


dG 

dv k 


2 n Vj G 


= 0 . 


4. Spectral Method Formulation 

Since we wish to consider a fully three-dimensional velocity space, we will set d = 3. When this 
three-dimensional velocity space (“v-space”) is coupled with a three-dimensional position space, 
the result is the full six-dimensional phase space of the single particle distribution function. In a 
three-dimensional relativistic velocity space, the domain is the interior of the sphere v 2 = 1 . so 
that a spherical polar coordinate system in velocity space is appropriate: 

(4.1) v = r[sin0(cos(px + sin(py) + cos0z] 


Here, x, y and z are orthononnal Cartesian unit vectors. In the velocity-space spherical 
coordinate system (v. 0. (p). equation (3.1 1) takes the form: 


dG 

dG 1 

1 v 2 f j(E + vxB) aG E-v 

dG / •>)!] 

— H V ■ — l 

v +2ml V“ Lr 

dt 

dx ' 
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The density integral (3.10) and corresponding current integral become (here, the element of solid 
angle is dL i = sin0z/0z/tp) 


(4.3) p = 



Gd\ = 1 - Jv 2 (l-v 2 ) 


GdvdQ. 
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(4.4) 



xGdx =- Jvv 2 (l-v") 


G dv dil . 


Our ultimate goal is to simulate numerically the evolution of the distribution function. To this 
end, a spectral method algorithm will be developed in which dependent variables are expanded 
in terms of known functions. Since the domain in v-space is a three-dimensional sphere, an 
appropriate choice of expansion functions for the angular dependence of the G is the well-known 
spherical harmonics T„„,(e.<P) [28] (the P„ m are the associated Legendre functions): 


(4.5) K ww (0,cp) = (-1) 


m 


j2n + l (n-mV 
4 k ( n + m ) ! 


cosey' ffMp . 


A choice of radial expansion functions is complicated by the relativistic factor ( 1 -v ) 1 in the 
integrands of (4.3) and (4.4). This difficulty is eliminated by choosing n = 5/2, which removes 
the relativistic factor from (4.3) and (4.4). although not from (4.2). [If n = 2 is chosen, the 
integrands in (4.3) and (4.4) appear to gain a weight function appropriate to Chebyshev 
polynomials. However, the domain of v is 0 < v < 1 , while the Chebyshev polynomials T n {z) are 
orthogonal over the interval -1 < z < 1, and a transformation z = 2v-l produces a weight function 
no longer in accord with the orthogonality properties of the 


Using n = 5/2 reduces (4.3) and (4.4) to simpler forms: 
1 ji 2 ji . 

(4.6) p = 1 - ]dv J dQ \ dty G v sin 0 

v =0 0=0 ( p =0 


1 n 2n . 

(4.7) j = - ) dv j dQ j Jcp G v v sin 0 . 

v =0 0=0 9=0 

The choice of radial expansion function is now straightforward. We will use a set of polynomials 
that are orthogonal over the interval 0 < v < 1 , with a weighting function of v‘. These are the 
shifted Jacobi polynomials of index (0,2): P„ 02 ( 2v-l), which are described in detail elsewhere 
[29]. Here, the necessary results are 
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The symbol P' n = dP n / ds denotes the first derivative of the /7 th order polynomial Pn(s). 
Using (4.9), let us define the normalized polynomial R„(v): 

a) R„(v) = v/2rt + 3 P®' 2 (2v - 1) 

(4.11) 

1 o 

b) \v~R n {v)R m (v)dv = &„ . 

0 


The distribution function G then has the expansion: 

00 00 f ] 

(4.12) G(v,0,(p) = I I ZY4m»^(v)K ;wi (e.q>)- 

k =0 n—{)ni~—n 


Here, it should be noted that the y* nw = yi are functions of position x and time t. 
The coefficients y knm can be determined by the inverse operation (■■■) knm ■ 

w.i3) y<™, ={c) ( .,„„ =JG(v.e, ( p)R t (v)y*je.< f »v 2 c (v<«. 


If (4. 1 2) is placed into (4.6) and (4. 1 3) is used, the result is 
(4.14) p = 1 — Ycx )0 - 
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To evaluate (4.7), (4. 1 ) can be written as 


4tt / * * * ^ 

(4.15) v = J— v(y n e + + T] _i e_ + Eioe 0 


Here, the complex unit vectors are 

x + iy 


(4.16) e + = — 


4~2 ' 


„ x-iy 

e_ = — 7 =^-, e n =z. 


V5 ' 


Similarly, using (4.1), (4.12) and (4.13), allows the integral (4.7) to be determined: 


(4.17) j = Jk Re 


VT5 


7 n 1+7 on 
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It is clear that the choice of expansion functions leads to an easy evaluation of the charge density 
and current, once the expansion functions of the distribution function are known. When p and j 
are determined, Maxwell’s equations (2.1 1) to (2.14) can be used to find E and B. These, in 
turn, affect the evolution of the y through (4.1). The spectral form of the evolution equations 

will be described next. 


5. Spectral Form of the Vlasov Equation 

The Vlasov equation (4.2) is, in a three-dimensional velocity space spherical polar coordinate 
system. 


(5.1) 




v 1 1 — v 


Xfac 
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-5vG 
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Here D is the angular derivative part of the velocity space gradient operator. 


(5.2) 


In the 

Now, 


= v — + -D, D = + 

9v 3v v 36 sin 0 dip 

above, v,0,cp are the unit vectors in the v, 0. <p directions, respectively. 

we subject equation (5.1) to the inverse operation (•")^ 7im defined in (4.12). 


(5.3) 
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where 


< 5 - 4 > <,!,H v c> 
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At this point, we use the properties [28] of the angular momentum operator, L: 

(5.8) L = -/vxD = -l=(L + e_ -L_e + ) + L-e 0 , 

V 2 


where 


(5.9) L+ = ±e ±/<p 
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— ±/cote^- 
30 3<p 


L. = -i 
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We will also need 


(5.10) D = -i v x L =-)=rcos0(L + e_ + L_e , )-e 0 sin0 — 

4i + + ° 30 


(5.11) v = -|=sin 0(e ,<p e _ -e ,tp e + )+ cos0e o . 
V 2 


The operators v , D, and L act on the spherical harmonics T„„ ; (0,<p) to produce linear 
combinations of other spherical harmonics with possibly raised and/or lowered indices n and/or 
m. The coefficients of these linear expansions can be made more succinct using 


10 



c 


m _ 
n ~ 


(// - m){n + m) 
(2n-\)(2n + \) 


( 5 .i 2 ) s;;' = 


(n + m )(/? + m — 1) 

]j (2//-D(2/?+ir 


T” 1 = *J(n-m)(n + m + \) 


Using the tabulated properties [28] of the JWG.tp), along with the above definitions, produces, 
after a lengthy derivation, the following representations of (5.4), (5.5), (5.6) and (5.7): 
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The terms u| a j above are 


(5. 1 6) w'j, 0) - = j v 3 Rj dv 


(5. 1 7) = J v 2 (l - v 2 f 2 (r) - 5Rj K ^ 
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We now have a set of partial differential equations (5.3) in position space (x-space) that is 
nonlinearly coupled to the Maxwell equations. The velocity dependent part of the Vlasov 
equation has thus been transformed into a set of v-space spectral equations, i.e., equations for the 
spectra Y*,„„(x,/). This v-space spectral formulation is independent of the numerical method that 
is used in x-space. In x-space. we are free to choose whatever geometry, coordinate system and 
numerical method seems most appropriate and realizable. In fully ionized plasma problems, 
there are no distinct bounding surfaces since these generally preclude the existence of a fully 
ionized state. We must, of course, make some reasonable assumptions about the behavior of the 
plasma at the numerical boundary. We could assume that the number density falls to zero at 
such a boundary, for example, or we could assume that density is periodic at the boundary of the 
numerical grid that represents x-space. This choice then narrows the choices of numerical 
method to use in x-space. Here, for illustrative purposes and because it is applicable to so-called 
homogeneous plasmas, we will choose periodic boundary conditions. This leads to a 
straightforward spectral formulation in x-space. Thus, our next step is to consider solution of the 
Maxwell equations in the case of periodic boundary conditions. 


6. Determination of the Electromagnetic Field 

The full set of Maxwell’s equations, (2. 1 1 ) to (2. 14), can be written in terms of the potential 
(p and the potential vectors A and C as follows [30]: 
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( 6 . 1 ) 



(6.2) I^ = V 2 A + 4;rj 
c dt 

(6.3) V 2 <p = -4;r/9-V-C 

(6.4) -L^ + VA = 0. 

c 2 dt 

Equation (6.4) is the Lorentz condition [31]. The relations (6.1-4) are a set of coupled, first- 
order-in-time, partial differential equations for the electromagnetic potentials, in a form suitable 
for numerical integration. The associated electric and magnetic fields are 


(6.5) B = Vx A 


(6.6) E = —V <p — C. 

In the case of periodic boundary conditions, we will expand electromagnetic fields and 
potentials, as well as the y klw » when necessary, using fast Fourier transforms: 

a) A(x,0 = (27tr^IA(k,r)c ikx 

k 

(6.7) 

b) A(k,/) = (27t) _ ^IA(x,t)c' lk x . 

X 

There is no time-frequency expansion here because the coupling between the electromagnetic 
field and matter is generally nonlinear, which generally precludes a simple dispersion relation 
between wave frequency w and wave vector k. In terms of Fourier coefficients (omitting t in the 
argument, for brevity), equations (6.1) to (6.6) take the form. 


( 6 . 8 ) 


1 dA(k) 

c dt 


= C(k) 


(6.9) - J ^ (k - = -A- 2 A(k) + 4^-j(k) 

c dt 
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(6.10) #>(k) = k~ 2 [47rp(k ) + ik -C(k )] 


(6.11) -L— + ik-A(k) = 0. 

c~ dt 


The sources in (6.9) and (6.1 1) are the Fourier transforms, using (6.7), of (4.14) and (4.17). 
(Also, note that the source coefficients corresponding to k = 0 are zero.) 

Using (6.7) and (6.10), equations (6.5) and (6.6) become 


(6.12) B(k) = ik x A(k) 


(6.13) 


E(k) = —47t i 


-^rP( k )- 
k 2 


kk 

~ 1 ? 


\ 

C(k). 

/ 


Here, I is the unit dyadic. Note that ^(k)does not appear in these expression for B(k)and 
E(k) . The scalar potential is thus only an auxiliary function, while A(k) and C(k) are seen to 
be more fundamental in classical electrodynamics. Furthermore, it is only those parts of A(k) 
and C(k) transverse to k that are essential to determining B(k)and E(k). (However, in 
numerical simulations of quantum electrodynamic processes, ^(k)must also be calculated 
[32,33].) 


7. The Complete Spectral Algorithm 

A spectral method algorithm for the Vlasov-Maxwell system, as developed here, consists of a 
part in velocity space and a part in position space. The basic v-space spectral equation is (5.3) 
and its solution produces the coefficients y („„„ which are time-dependent functions in x-space. 

We use a small subset of these coefficients to directly determine the electromagnetic sources, 
through (4.14) and (4.17). In x-space, Fourier transforms of these sources then allow the time- 
evolution of the electric and magnetic fields to be calculated through (6.8), (6.9), (6.12), and 

(6.13). These electromagnetic fields then enter the right-hand side of (5.3), completing the cycle. 
This is a coupled set of nonlinear equations, and a complete description of the algorithm 
developed here requires only a few additional details. 
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Equation (5.3) determines y knm and can be written as 


(7.D = 

dt 

where 

(7.2) r taI =-VN ™»,. n™, =E'kl + 

To determine the r km „ we use the Fourier transform to k-space (6.20b), to find N^,, and then 
form 

(7.3) r knm --tk-N A7IW) , 

after which we transform back to x-space with (6.20a) to get T klwl . Similarly, using (6.20a) to 
find E(x) and B(x) from E(k) and B(k) , allows us to form point-wise on the numerical x- 
space grid. The transform method was pioneered many years ago [6] and has proven very useful 
in the solution of nonlinear problems involving many modes. 

Equation (7.1) is now defined on all x-space grid points and can be used to advance the y knm 
forward in time with a suitable time-integration scheme (such as a Runge-Kutta [34] or an 
Adams-Bashforth [35] procedure). Since r knm consists of only six coefficients y knm , it may be 
efficient and more accurate to treat this term implicitly when solving (7.1). Aliasing (i.e., higher- 
order coefficients affecting the determination of lower-order coefficients) can occur in 
determining the H\„„„ but this can be eliminated, if desired [36], However, aliasing is often not 
critical, and run time can be reduced by neglecting aliasing effects (such methods are termed 
pseudospectral [37]). Lastly, one can introduce external parts to E and B in (7.2), in addition to 
those parts determined by (6. 1 2) and (6. 1 3), and these can serve as a source of external forcing, 
if appropriate to a given problem. 


8. Discussion 

The spectral method algorithm developed here provides a means for solving the coupled, 
nonlinear Vlasov-Maxwell system of equations in the full six-dimensional case. Numerical 
implementation requires the use of a discrete grid of points, and if we assume, for the purpose of 
illustration, that each dimension is given N points, then there are a total of N 6 grid points. \fN = 
32, then /V 6 - 10 y , and if we assume 8 bytes per word in a computer, along with a factor of 25 to 
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50 for all the arrays needed, then core memory requirements are in the range of 1 to 2 terabytes 
( 10 1- bytes), which is currently available on the largest existing supercomputers. The number of 
grid points allotted to each dimension in the six-dimensional x-v-space can also vary, to increase 
or decrease resolution, as appropriate. 

However, remember that, in a spectral method, grid points are sampling points for an underlying 
continuous representation. The number of grid points is thus essentially equivalent to the 
number of coefficients (the ‘spectra’) kept after truncating a formally infinite function expansion. 
To prevent a ‘pileup’ at the higher-order coefficients and to ameliorate the effects of 
‘fi lamentation’ [10,12,17], it may prove necessary to investigate the use of a collision term in the 
Vlasov equation, as done previously [3,16]. (These possibilities, and others unthought of, remain 
as potential challenges for future development.) 

In spherical polar coordinates in v-space, the number of points in the radial, or v, direction sets 

the energy resolution of the code. This resolution may be defined by the relativistic factor y„ = 

-> -V2 

( l-v„“) corresponding to the largest positive zero v„ of R„(v) (all the zeros are used in the 
Gaussian quadrature formula [38] for integrals involving R„, in which case the function to be 
evaluated is effectively known more or less precisely only at those zeros and interpolated in 
between [39]). Using known results [29], we can estimate y„ ~ &n/3n. For n = 32, this gives the 
resolution of electron energies an upper bound of around 14 MeV, a value which can be 
increased by making n larger (c.g., by decreasing the number of grid points allotted to other 
coordinates). Thus, moderate resolution six-dimensional numerical simulations of relativistic 
electron distributions are possible on available computers. 


9. Conclusion 

The algorithm presented here provides a straightforward procedure with which to begin using 
supercomputer resources. The next steps are to evaluate certain integrals mentioned in this 
paper, to begin the numerical implementation of the algorithm, to gain access to the necessary 
computers, and to begin running numerical simulations. The probable outcome of proceeding is 
the creation of a numerical tool capable of studying the full six-dimensional dynamic evolution 
of charged-particle distributions, starting from various initial conditions and taking into account 
the influence of external fields. The benefits to be gained are a greater understanding of basic 
plasma processes and a new capability to predict radiation levels and effects in astrophysical, 
spaceflight, and engineering systems. 
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